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The  subroutine  for  the  calculation  of  intensity  deviation  in  Lowtrsn  V 


I.  Introduction 

In  this  report,  we  present  theoretical  formulas  and  a  coded  subroutine 
for  the  calculation  of  intensity  deviation  which  shows  the  high  bound  and  low 
bound  of  the  atmospheric  transmittance  due  to  the  turbulence  effect.  In  the 
subroutine,  we  calculate  the  transmittance  deviation  for  point  receivers  as 
well  as  for  finite  aperture  receivers  which  exhibit  the  aperture  averaging 
effect. 

Depending  on  the  type  of  propagation  paths,  the  calculation  in  the 
subroutine  has  been  divided  into  three  parts -horizqptal ,  upward  and  down¬ 
ward  paths,  similarly  to  the  transmittance  calculation  in  Lowtran  V  program! 
The  different  examples  will  be  shown  and  compared  in  the  following  section. 


II.  Plane  wave  intensity  deviation  for  point  receivers 

Consider  a  plane  wave  U  propagating  through  the  turbulent  medium  repre¬ 
sented  by 

U  -  eX+  iS  (2.1) 


where  x  is  a  real  value  that  represents  the  log-amplitude  and  S  is  the  imagin¬ 
ary  part  that  represents  the  phase.  Assuming  the  Gaussian  probability  distri¬ 
bution  for  x,  the  average  intensity  and  variance  of  intensity  can  then  be 
stated  as 


(2.2) 


<12>*(U'U*U'U*)  -  e4<*'  + 


(2.3) 


2  2 

where  (  )  denotes  an  ensemble  average  and  c  2  ((X  *  (x))  ^sing  Eq.  (2.2) 
and  (2.3),  the  normalized  variance  of  the  intensity  fluctuations  is  given  by 

c*  »  ilh  -  W2.  .  .(  (2.4) 

1  <I>2 

2 

The  variance  of  log -amplitude,  ,  has  been  found  by  Rytov’s  Method.  However, 

it  is  only  valid  for  weak  turbulence  when  applied  in  Eq.  (2.4).  Experimental 

data  indicates  that  cr^  is  saturated  toward  the  value  of  unity.  In  recent 

years  theoretical  work  to  proof  that  the  variance  of  intensity  saturates 

2 

to  a  constant  of  unity  was  performed.  Avoiding  complex  mathematics  and  hoping 
to  get  a  model  which  is  sufficiently  accurate  under  weak  and  strong  turbulence 
conditions,  we  relate  the  variance  of  intensity  and  log-amplitude  by 

a  «  1  -  e'2a-c  (2.5) 


For  small  values  of  o  ,a  «■  2  a  which  agrees  with  equation  (2.4).  For  large 

X  1  X 

a  ,  ffj  »  1,  which  agrees  with  the  saturation  condition.  After  we  determine 
X 

Oj.,  the  upper  bound  and  lower  bound  of  transmittance  can  be  calculated  by 


T 


l 


»  T(1  +  Jj.) 


(2.6) 


using  Ref.  (2),  the  variance  of  log-amplitude  as  found  by  Rytov’s  Method  is 
given  by 


CT2  *  .563k7/6  rL  dTl  C2  (n  )  (L  -  Tl>5/6  (2.7) 

X  J  o  n 

2 

where  k  is  the  wavenumber,  L  is  the  path  length  and  is  the  structure 

2 

constant  of  refractive  index.  When  Cn  is  constant  along  the  path,  Eq.  (2.7) 
can  be  rewritten  as 


,31  C2  (h)k7/6LU/6 
n 


(2.8) 


2  5 

The  structure  constant  Cn  has  been  measured  and  modeled  by  Hufnagel,  et.  al. 


we  have  modified  it  to  fit  in  Lowtran  as 


Cn  <h> 


4.2  x  10* 14  h'273  exp  (-h/ 320) 
3.77  x  10‘15 
0 


(h  >  10m) 
(h<  10m) 

(h  >  100km) 


(2.9) 


where  h  Is  altitude  and  is  in  units  of  meters. 

For  horizontal  path,  h  is  constant,  hence  we  use  Eq.  (2.8).  For  downward 
(and  upward)  path.wemust  use  Eq.  (2.7)  in  which  the  integral  shows  an  in¬ 
tegration  from  transmitter  to  receiver.  From  the  weight  function  (L-1))5^ 
in  Eq.  (2.7)  we  know  the  turbulence  around  transmitter  has  more  effect  than 
the  turbulence  near  the  receiver.  Hence  we  predict  that  a  downward  path  has 
larger  variance  of  intensity  than  an  upward  path  for  the  same  length  of 
path. 


For  the  program,  we  rewrite  Eq .  (2.7)  in  summation  form  as 

2  ...  7/6  _  ..  2 


a  -  .56k?/6  E  1  (h  )  (L-L  )5/6  — - - - 

X  i  j  iJ  lj  hi  ’  hi-l 


Ah 


iJ 


(2.10) 


hj^<25km  A  *  20m 

25km  <hij<50km  A  h^  -  100m 

50km  <  h^j  <  100km  A  h^^  =*  400m 

where  hjj  is  the  altitude  corresponding  to  the  calculated  point  of  the  path, 
"i"  is  the  layer  index,  "j”  is  the  sub  index  of  each  layer,  L  is  the  total 
path  length,  Lj.  is  the  path  distance  from  transmitter  to  the  point  calculated 
and  ALj  is  the  path  length  for  each  layer  passed.  The  choice  of  Ah  intervals 
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was  made  to  allow  better  height  resolution  in  future  speci fications  of  C“(h). 

n 

In  this  calculation,  we  assume  that  the  path  in  each  layer  is  straight. 
Refraction  occurs  only  at  the  boundaries.  The  refraction  calculation  is 
executed  by  the  original  Lowtran  program. 


Figs.  (2.1),  (2.2)  and  (2.3)  show  the  transmittance  predicted  by  Lowtran  V 
with  the  deviation  calculated  by  the  added  subroutine  for  horizontal,  downward 
and  upward  path,  respectively.  These  calculations  are  for  point  receivers  and 
a  5km  path  length.  We  find  that  the  deviation  for  downward  path  is  larger  than 
that  for  upward.  This  is  due  to  the  stronger  turbulence  at  the  lower  altitude 
and  the  effect  that  turbulence  near  the  transmitter  is  dominant. 

III.  Intensity  deviation  for  finite  aperture  receivers 

In  Sec.  H»the  intensity  fluctuations  are  assumed  to  be  measured  by  a  point 
receiver.  However,  in  the  real  world  the  receiver  has  a  finite  aperture.  If 
the  diameter  of  the  objective  is  much  larger  than  the  amplitude  correlation  dis¬ 
tance  p,  the  objective  will  contain  wave  front  sections  with  fluctuations  of 
opposite  sign,  so  that  the  overall  light  flux  through  a  larger  objective  will 
fluctuate  relatively  weakly  compared  to  the  flux  through  a  small  (compared  to 
p)  objective. 

Consider  the  intensity  in  the  receiver  plane  to  be  I  (L,  £  ) .  The  total 

flux  P  through  the  objective  is  then  P  *  J  I(L,  £)d£,  where  £.  is  the  transverse 

T 

coordinate  and  Y.  is  the  aperture  area  of  the  objective.  The  fluctuations  in  P, 
defined  asP^P-^P),  are  expressed  in  the  form  P'=  jjj  l'  (L,£)  d£,  where  I1 
I  -  (  I  For  the  mean  square  fluctuations  of  power  we  have 

<p/?>  s  <  I#  <Ri>  T'  <E.Z1>  dEl  dE-2 

Y  Z 

00 

=  J"J  B1  (£l  ■  •E2)  F(s1)  F(£2}  d£i  d£2  (3a) 

-00 

where  the  intensity  correlation  function  s  (i'  (g.^)  I'Cg^)'1’  and 

F(P)  is  a  function  which  is  zero  outside  the  aperture  and  1  on  its  surface. 


range  5  km  for  a  point  receiver  (D  =  0). 


receiver. 
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Changing  variables  and  assuming  the  receiver  aperture  to  be  a  circle  with  radius 
R,  Eq.  (3.1)  can  then  be  written  as: 


2R 


(P*  )  ■  2rr  pdp 


Bt(p) 


K(p) 


(3.2) 


where 


K(o)  H  J  r  F(£l)  F(£i-  o)  d2£l 


2R2  [cos  -1  (£)  *  -^1*  P<2R 


(3.3) 


0, 


>2R 


The  normalized  fluctuation  of  power  is  defined  as 


Q(R) 


,  p'2  ■,  .  _4_  ,2R 

<  p>2  m2  0  <i>2 


[cos 


-1 


(_C_) 

2R; 


(3.4) 


to  show  the  averaging  action,  we  define  another  parameter  G(R)  =  which 
is  a  ratio  that  compares  the  fluctuations  of  power  of  a  finite  aperture  and  a 
point  aperture: 


G(R) 


(P) 


[  cos 


]  pdp 


(3.5) 


where  bj(p)  is  the  normalized  correlation  coefficient  of  the  intensity  fluct¬ 
uations 


bj(p)  ® 


bt(p) 

—■ T  4B.,(p) 

HZ 


Br(0) 

<T>2 


-  i 


(3.6) 


4<j 


-  1 


B  (p)  =  <X  (E.1)x  '  (22)  > 


and 


(3.7) 
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2 

Because  the  structure  constant  depends  on  the  altitude  h,  the  calculations 
for  horizontal  and  downward  or  upward  paths  are  different.  We  separate  the 
two  cases  as  follows: 


(i)  Horizontal  path 

By  using  Rytov's  Method  and  the  Kolmogorov  spectrum  (Refs.  2,6),  B (p)  can 

2 

be  found  with  the  condition  that  -t  <  <  \L  as 

o 


and 


1  -  12.3 


V»>  ’  VB><5x 

2 

_ 2 _ 


(3.8) 


au5\m 


P<  <tr 


bx<p)  -  ^  1  -  2.36(-a-L---)5/6  +  1.71 

kp2  '7/6 
-  .0242  (^ 


2  2 

.024  (&-  --)  ^  lQ«  p  «  /Tl 


,  J~\ L<< 

(3.9) 


where  4^  is  the  small  scale  of  turbulence  and  is  assumed  to  be  3mm  in  the  Low- 

tran  code,  and  ct  is  the  same  as  in  Eq.  (2.8).  After  substituting  Eq.  (3.8) 

X 

(2.8)  and  (3.6)  into  (3.5),  we  can  find  the  receiver  aperture  averaging  effect 
for  horizontal  paths. 

(ii)  Slant  path 

Using  Rytov  Method,  Kolmogorov  spectrum  and  locally  homogeneue  medium,  the 

correlation  function  of  log-amplitude  can  be  found  from  Ref.  (  2  )  under  the 
2 

condition  L\  >>  : 

B  (L.p)  -  .033  tt2(-  F(  -  )  k2*  ;L  C2  (T))  dT)  •FG(p.Tl)  (3.10) 

X  D  o  “ 


where 

FG(o,D) 


\2  k  1  lV  6, * 


- E - ) 

* 


10 


1 5/6 
m 


F 

11 


('  6  ’1> 


2  2 

x.  p 

m 


(3.11) 


k  .-5*2L 

®  i 

u 


(3.  2) 


^F^(a,b,x)  is  the  degenerate  hypergeometric  function  .  For  calculat  .ng 
on  the  computer,  we  approximate  the  function  ^F^  as 


1  -  0.8333X  -  0.0347X2  -  0.0045X3  |x|<6.5 


!F1  <-f  1>X) 


1 .0627 (  -  X) 


5/6 


(3.13) 

Ixl  2  6.5  and  (Re  X  <  0) 


This  approximation  has  been  checked  arid  the  total  error  is  under  10%.  Further¬ 
more,  we  assume  -  ^  .  .1  ..Tl .  ;>  > — i—  which  means  that  we  neglect  very  small  parts 

m 

of  turbulence  near  the  receiver.  We  then  substitute  Eq.  (3.13)  into  (3.11), 
and  get  the  different  forms  for  the  different  conditions: 

2  2 


D5/6(.259  +  .8U5  -£ 


FG(71,pW  -  1.065(p2/4)5/6 


.5/6, 


4D 


2 

4D 


4D 


1r>2- 


<45- >6-5’ 

4A 

>6.5  ) 

2  - 

0043  (-^)J  ) 

<-4D-<6*5’ 

2 

p 

•%  z  c\ 

4A 

b •  j  J 

2  - 

.0043  (-^)J) 

(3.14) 

2  2  2 
-  (A)  5/6  (1  +  .8333  -  .0347(-^~)2  +  .0045(-J^-)3) 


( 4F<6-5-  ir 


<  6.5  ) 


#Note  that  an  equation  for  a  similar  situation  given  by  Ishimaru  (A.  Ishimaru, 
"Wave  Propagation  and  Scattering  in  Random  Media",  Academic  Press  1978,  Vol.  2, 
Eq.  17-112)  is  in  error. 
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where 


,  A  s 


(3.15) 


Putting  Eq.  (3.13)  into  (3.10)  and  changing  the  integration  into  summation  we 


2  2 

V)*2a7k  l  *  C„  (l‘lJ>re«-tJ> 


(3.16) 


where  k  is  the  wavenumber,  L  is  the  path  length,  "i"  is  the  layer  index,  "j" 

2 

is  the  subindex  in  each  layer,  C  (h,  .  )  is  the  structure  constant  at  the  height 

n  lj 

j  ’  ^ij  C^e  a^t^Cu<^e  °f  Che  point  calculated,  and  is  the  path  length 
in  each  layer.  Similarly,  putting  Eq .  (2.10), (3.6) ,  (3.16)  into  (3.5)  and 
changing  it  into  summation  we  get 


G(R)  »  I  bJ.(2Ryi)  [cos‘1(yi)  -  y^J  1  -y*  ]  y t  *.01 


(3.17) 


where  "R"  is  the  radius  of  the  receiver  aperture,  i  is  the  summation  index, 
and  y  =  O.Oli. 

The  intensity  deviation  a^(R)  for  a  finite  aperture  receiver  can  now  be 
obtained  by  multiplying  G(R)  and  o  (Eq.  2.5) 


a  (R)  -  o  *VG(R) 
P  1 


(3.18) 


2  -16  -?/3 

For  a  condition  of  moderate  turbulence  C  =*2-  10  m  ”  and  horizontal  pro~ 

n  r 

pagation  with  path  length,  L  ■  5km,  the  aperture  averaging  coefficient  G(R) 

is  shown  as  solid  line  in  Fig.  (3.1).  It  is  similar  to  the  dashed  line  2  in 

2 

Fig.  (3.1)  predicted  by  Tatarskii  .  The  dashed  line  1  in  Fig.  (3.1)  shows 
the  aperture  averaging  coefficient  G(R)  for  downward  path  with  altitude  from 
200m  to  4km.  Comparing  horizontal  and  downward  cases,  we  see  that  aperture 
averaging  has  more  effect  for  the  horizontal  case.  This  is  because  the 
total  turbulence  in  the  case  of  horizontal  path  is  stronger  than  the  downward 
case  and  the  coherence  length  of  field  of  latter  is  longer  than  the  former. 
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Figs.  (3.2),  (3.3)  and  (3.4)  show  the  atmospheric  transmittance  predicted 
by  the  modified  Lowtran  V  with  diameter  of  the  aperture  being  30cm,  for  hori¬ 
zontal,  downward  and  upward  path,  respectively. 

IV.  Subroutine  of  intensity  deviation  for  Lowtran  V 

This  subroutine  is  for  the  calculation  of  intensity  deviation,  due  to 
turbulence,  which  can  be  used  to  define  the  upper  bound  and  the  lower  bound 
of  plane  wave  transmittance  for  the  point  receiver  case  and  the  finite  aperture 
receiver  case.  The  subroutine  for  each  calculation  of  transmittance  and  for 
each  frequency  is  called  by  subroutine  TRANS,  one  of  the  subroutines  of  the 
main  program.  According  to  horizontal  path,  downward  path,  upward  path, 
we  divide  the  subroutine  into  three  parts.  The  attached  flow  chart  shows  the 
main  points  of  the  subroutine  (See  Appendix  A).  We  transfer  the  data  of  path 
altitudes  HI  and  H2,  path  length  L,  wavenumber  k,  the  path  length  in  each  lay¬ 
er  DS1  ,  DS2,  the  height  in  each  layer  XW1,  XW2,  and  the  diameter  of  the  rece¬ 
iver^  and  other  necessary  data  to  the  subroutine  VRANI  from  the  main  program 
or  from  the  subroutines  of  the  main  program.  The  symbols  and  definitions  of 
variables  used  in  the  subroutine  are  listed  in  Appendix  B.  The  program  list¬ 
ing  is  in  Appendix  C. 

V.  Limitations  and  comments 

The  subroutine  for  intensity  deviation  is  theoretically  formulated  by 
Rytov’s  He  : hod.  Even  after  our  modification,  Eq.  (2.5)  is  only  valid  for  weak 
and  extremely  strong  turbulence.  For  moderate  turbulence  the  theory  under 
estimates  the  fluctuations  of  intensity  as  compared  with  experimental  data. 
Second,  in  this  program  we  only  consider  the  plane  wave  case,  hence,  it  will 


#The  diameter  of  the  receiver  must  be  read  in  from  the  sixth  variable  of  the 
second  control  card  and  the  format  is  F10.3  and  in  units  of  meters. 
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overestimate  the  fluctuation  In  weak  turbulence  for  a  real  source  such  as 
a  beam-wave  or  spherical  wave. 

For  the  turbulence  itself,  we  use  the  model  given  in  Ref.  5,  which 
seems  to  be  too  simple  to  give  satisfactory  result.  Hopefully,  we  can  find 
a  more  proper  form  of  turbulence  profile  in  the  future  and  update  it  in  the 
code.  In  fact,  turbulence  in  the  atmosphere  should  be  varying  with  tempera¬ 
ture,  wind  speed  and  constituents  of  the  atmosphere.  It  is  suggested  that 
the  turbulence  model,  like  the  aerosol  models  or  atmospheric  models  in 
Lowtran,  should  be  constructed  by  different  models  based  on  different  areas 


and  seasons. 


AVERAGE  FACTOR  GfR) 


receiver  rrdius  [cm] 


RfiDIUS/FRESNEL  HONE  (R^T 1 


Fig.  3-1  Comparison  of  aperture  averaging  factors  G(R)  for 
different  situations.  Solid  line  is  predicted  by 
the  new  subroutine  for  a  horizontal  path  at  altitude 
h  =  400m.  Dashed  line  1  ( —  -  —  )  is  predicted  by 
the  new  subroutine  for  a  downward  path  at  altitude 
from  200m  to  4km.  Dashed  line  2  ( —  — )  is  pre¬ 
dicted  by  Tatarskii.  The  wavelength  is  X  -  1^  m 
and  path  length  is  L  *  5km. 


receiver  of  diameter  30  cm. 


L OUTRAN  V  PREDICTIONS  WITH  VARIANCES 
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finite  aperture  receiver  of  diameter  30  cm. 
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to  200  m  with  visual  range 
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APPENDIX  A 

Flow  Chart  of  Subroutine  VRANI 


APPENDIX  B 


VRANI  Symbols  and  Definitions 

AD  Diameter  of  receiver  aperture  in  meter  (m) . 

ANGLE  Initial  zenith  angle  in  degree 

BI  Covariance  of  intensity 

BL  Normalized  covariance  of  log-amplitude 

BX  Covariance  of  log-amplitude 

2 

CN2  C  -  structure  constant  of  turbulence 
n 

DD  the  ratio  of  the  distance  from  point  calculated  to 
receiver  over  total  path  length. 

DH  Height  interval  of  slant  path  integration 

DS  Distance  from  point  calculated  to  receiver 

DT  Same  as  DS,  especially  used  in  the  downward  long 
path  calculation 

DZW  Difference  of  height  between  two  near  layers 
FR  Fresnel  zone  in  meter  (m) 

GD  Aperture  averaging  coefficient 

HMIN  The  minimum  height  of  a  downward  path 

HW  Height  corresponding  to  the  point  calculated 
HI  Height  of  transmitter  (and  receiver  for  horizontal 
path) 

H2  Height  of  receiver 

IV  wavenumber  in  cm  ^ 

JMIN  the  layer  index  of  the  minimum  height  for  a  downward 
pith 
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RANGE  Path  length  in  kilometer  (km) 

VR  Variance  of  log-amplitude 

VRI  Intensity  deviation 

PVR  Power  variance 

DSW  Path  length  in  a  layer 

WH2  Height  of  receiver  in  meter  (m) 

WL  Wavelength  in  meter  (m) 

WK  2tt/\  ,  wavenumber  in  m  ^ 

WRANGE  Path  length  in  meter  (m) 

WV  Intensity  variance  without  approximation 

XWl  The  lowest  height  of  a  given  path  in  a  given  layer 

XW2  The  highest  height  of  a  given  path  in  a  given  layer 
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APPENDIX  C 

Listing  of  Subroutine  VRANI  Program 


5u9»0u7InE  vSanI(IV) 

*•  *  vR 

C  THIS  SUdROuTInE  13  TO  CALCULATE  THE  VARIANCE  OF  INTENSITY  avB 

C  OUE  TO  TURBULENCE  ANO  TH£  CALCULATED  STANOARO  UEVIATION  CAN  BE  AVR 

C  USED  TO  DEFINE  high  BOUND  And  LUn  SOUND  OF  TRanSMI  TTanCE  avR 

C  AvR 

COMMON  /CARO  1 /  HODEL, IhaJE, ITYPE, lEN, JP, IH.hi , H2, mj, ml. lE^ISS, RO  avr 
1  .TdOUNO, IStASN, IVULCN.vlS  AvR 

COMMON  / CARoa/  Hi.Hj, ANGLE. range, beta, MHlN, RE, au  avr 

COMMON  /CAROS/  vi ,  Vi,DV,AvM,COiCM,«aS]  ,£(151  ,CA,f>!  AVR 

COMMON  /CNTRL/  LENST.ANAA.M, IJ.Jl, ja.JMIN, JErTRA, IL, IKMAX.NLL.NR1  AVR 
I.IFINO.NL.IKLQ  AVR 

COMMON  /WANG/  *a,0St.(J«)  .DZ«(3«)  ,A*1  C3R)  .**2(34)  AVR 

Common  /vm an/  vri  4Vr 

0IMEN3I0N  «3(3a)  avr 

H*FH1«1000,0  avr 

«RANG£iRAnGE»1000,0  avr 

«h2«mJ« 1000. 3  av£ 

VRA0.W  A  VR 

Cn2»0,0  avr 

Pv  R*0 . 0  avr 

ML«0,«l/IV  avr 

FR « (HL»mRANGE1 »>0.S  AvR 

MR»IV.t00,.i,.Pt  avr 

VK«MI<«  4(7. /6.)  avr 

40*9 ,E»a/5, 9 10*»a ,  avr 

IFCITYPE.NE.l)  go  to  20  AVR 

C  AVR 

c  variance  calculation  for  horizontal  path  avr 

C  A  vR 

CN2«R,aE-lR*MM««  (-2,/3.)«e*P(-MM/5ZtS,0)  AVR 

IFCht.GT, 130,3)  CN2«n,l  avR 

IF  iMM.Lt.10,0)  CNd»B,7TF-i5  avR 

VW “W, 3 1 •CN2»HMANG£*R ( 1 1 ,/Rt )  A  V  R 

VRRVHaVK  AVR 

VHl.lf.E*PI-a..VR..0,S)  AVR 

IFIaO.lt. U.00VI  GO  TP  9 1  AvR 

OU  IB  I«1.U)0  avR 

r«0,0i«i  avr 

Oy»»0«y  avr 

IF CUY, GE, 0.0833  GO  TO  U  avr 

BL»1  .•1J.J»OYAA0.0/ CFRF«Ib,/3,) •0,0B3RR(t  ,/3,  J )  AvR 

GO  TO  IT  AVR 

u  YiRRFAnYARa.a/HRANGE  av» 

IFIOY.GE.FR)  GO  TO  12  avR 

BL»l,«*2,3a«Jl»»(5,/6,)»l,71»«I»0t,J2««XI«»2,0  AVR 

GO  TO  17  AVR 

12  BLRM0,02«2Rim/4,)R«I-T,/B,J  AVR 

17  Continue  avr 

B*»t)L«V«  avr 

BI«£*P(4,0«tf X)»l,  AVR 

RVR.PVR*8I. (ACOS(Y)-Y.(l..Y..i ,) a.0.5) AYV0.01  AVR 

18  CONTINUE  AVR 

PvHaPv w« 1 8 , /B I  AvR 

MVAt*P(4,0*V«)"l ,  A  V  B 

GO*PVR/.V  AVR 

VHI«VRI.GOaa0,B  avr 

GO  TO  91  AVR 

23  IF  C  ANGLE , ST, 9«  ,  0)  GO  r0  37  A  v  R 

i<» 
A  .  M 
A  V  R 


001 

002 

003 

004 

003 

008 

007 

006 

009 

eie 

011 

012 

013 

014 

015 

018 

017 

018 

019 

020 

021 

022 

023 

024 

025 

328 

027 

o28 

020 

u30 

03l 

032 

233 

034 

035 

03a 

037 

018 

039 

040 

041 

040 

0  “  3 

£44 

045 

046 

047 

048 

049 

050 
051 
052 
o5  3 
054 
055 
o5a 
o57 
058 
054 
^60 
06  1 
082 


C 

c 

: 


VAMlANCE  Calculation  Fy*  OjmnmaRO  HATH 
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■  >U  38  «»l,li.> 

1  VX 

003 

oi»o . « 

AvX 

004 

-5**y. .  0 

A  V  B 

000 

V»8.iU  *A 

AVB 

00  b 

0  v  *  a  0  •  r 

AVP 

307 

»E  »yt»«2 . /4 , 

AVP 

00b 

y  x  *  u ,  >j 

AVB 

300 

AVB 

0?0 

CU  34  1  > J l . J 2 

AVS 

071 

«S(I)*Obw(lj/uZ«(l) 

A  vS 

372 

IK1W1U)  .LE.2S.0J  GO  'o  J1 

AVB 

073 

IEUwi  (J)  .LE.5M.01  GO  TO  22 

AVB 

0  7  4 

nx»«o0 ,u 

AVP 

075 

GO  TO  23 

4  VP 

o7  0 

Un»20 , O 

A  VP 

377 

GU  TO  23 

AVB 

878 

UX*100,0 

AVB 

079 

•I*»t»W2(I)-<<0(I)  )  •10*0.0/ UH 

AVB 

083 

lAlMlli 

4  V  P 

081 

M«»XW) ( I ) «1 J0M.0 

AVB 

382 

JU  33  J*1  ,  I* 

AVB 

383 

n.any.OA 

4  VP 

384 

IE  [hw.r,t.O<»b20.0.»NU.  1  TYPE. £0,3)  XG2(1)  •100003,0 

4  V  P 

085 

l  r (H«,G£,OObO0,O)  GO  TO  34 

A  V  B 

380 

4  VP 

087 

vK»vP^^t5o*CN^*(k>iKAT>*Gt«0S)**(5t 

/b.)  ■OH*«SU) 

4  V  P 

088 

uij>  4N(,fc»uS)  /  •»* 

AVP 

080 

•?  » ►<  /  ij  0 

4  V  P 

800 

IE  <00. L r .-0)  GO  TO  27 

4  VP 

301 

IE (S.GE .h.S.AOO.O.Gt ,0,5)  GO  TO 

27 

4  VP 

002 

IE tx.LT.b.S.AM/.O.LT.b.S)  GO  TO 

28 

4  V  P 

003 

Gar>o»«(S./(,,)  •  (u,2SO*3,8U5*O*0. 

009«B«*2,»M,3043aBa«3,)-l,0b3«RE««(AVR 

004 

•S./o. J 

AVB 

005 

IE  Cg.lE.0,3)  G*0. 

A  V  P 

000 

go  re  j  i 

AVB 

007 

G«0.y 

4  VP 

008 

Go  TO  31 

A  VP 

009 

2  a  G*0„..  (5./0  ,  1  •  (0.2S9O  .  1  05»x.»0  ,0P9*.i  »«2.«0 ,  je  l  Ja0.«3.  J  (5  ,/0, )  AVP 


•  •ll,*«>,'>33Ja‘J*0,i.!347»ua»2.»O,0W4S«Q»»3.)  »  VP 

It-tli.LE.i3.il)  G»8.3  *VH 

SI  sx»bX*<.*02*Q*»*5U)  AvB 

0i»US*l>x««S(I)  ‘Its 

33  DUMlNUt  ‘VP 

m.»«k*MjE*Di  */B 

IE (UU.LE.0.)  oSaOb-OH«>o(I)  »V» 

vx«vP.*<,5o«CNi»(.«ff*'-<Gt-05)  •»<5,/O,)«CXl«2Cn»!000,-X'*)  »»S(I)  AVP 

8j(«(j*»CN0»r.«c*Le  1 1 ) « i  ?!*u,  *"«)  ««s  ( i )  *vb 

03«US*(A*2U)  *1000.0-mg)  (I)  **B 

3V  CONTINUE  AVP 

IE  (Ai).L  r. 0,001)  GO  TO  3b  *0B 

6X«a««2.  AVP 

3l«E*P(4, •«»)-!  ,  »»B 

Py  XiPvB.Xl*  (ACOS(T)»v»  1 1  ,-V»«2.)  »»W,5)  ‘ToB.at  *vR 

33  CONTINUE  ‘TB 

PyBaPvu*lb,/PI  A  V® 

3b  VHtiu.y.  AvB 

vxt»l.-E*P(-a.*yP«»«).3)  *VB 

IE  UO.LT. 0,001)  GO  TO  01  AVB 

*v  *fc lb (« . »vx ) •) ,  *  VP 

Gu*pvx /  »v  *  V* 

y « I avx I «GU*«L . 5  *vR 

GJ  TO  01  A*B 


?i3G 

>9101 

Jtiaa 

.1133 

3134 

3133 

313b 

0137 

3130 

0100 

3113 

3111 

3112 

3113 

3114 
*1 18 
311b 
3117 
31  18 
3110 

3123 

3121 

3122 
?  1 23 

3124 


23 


37 

L  L)  N  f  l\ut 

A  V* 

3i  as 

A  V<4 

atao 

VintANCf.  CaLuM.*  <1  J*<  F(J4  uXaAAU  ><»i« 

Al-i 

•naT 

A  id 

nao 

•1  u  bi  **»  1 1  111 0 

A  V  R 

1104 

lljld.c 

A  VR 

31  3d 

0  T  »w  ,  / 

AH 

3131 

OX  *0.  u 

A  VR 

?iia 

a  vR 

3133 

or«xo*r 

AvR 

3134 

9F«uT»«a,/a, 

A  V  R 

3135 

VH*3,d 

A  VR 

31  3b 

Q*«F/.u 

AvR 

3i  37 

L1«J1 

A  vR 

3138 

JU  bd  t*l  <•'*(. 

AvR 

3139 

01*1.1-1 

A  y  R 

31*3 

«Sa.U*0S»(kl)/'-'2«(LI) 

AVR 

3141 

if  <x»i  an  ,l£.<js,9)  oo  rc  38 

AVR 

3i«a 

if (*»t an .ue.sm.k-)  go  to  39 

AVR 

3143 

l)n*i»l'3,3 

AvR 

3144 

Go  TO  AU 

A  vR 

3  1  “3 

s 8 

0H*dt'  ,il« 

AVR 

314b 

Gu  10  4<i 

A  yR 

3147 

39 

0««  1  u«j .  oil 

AVR 

314b 

■«l**(Xwl(l.lj-5*atl.l31*lB*ia,d/0" 

AvR 

3144 

IKM, IK 

AVR 

3130 

( L  1 J  «13uH,J 

AVR 

3 1  S  1 

0u  57  J*t ,  t> 

AVR 

3153 

hnAHX-Om 

AvR 

■)  1  5  3 

CNam.it-i  *•  *.{-«,/  3,  )»txp(-H«/3aw.0) 

AVR 

Jt  5u 

v»»v&*a.5b*CNi* (««w;fc-os)  ••(%,/fe.)*OH«MscLi) 

AVR 

3155 

Du*  l-«iA(ifu3)  /  a  H 

AVR 

3 1 5b 

3  *RF  /t>D 

avR 

315  7 

IF (uO,LT.aOJ  GO  to  si 

AvR 

0  1  5b 

IF  GO  TO  5 1 

A  VR 

3159 

IF (x.lT.fi.5,»NU,u.LT,b,5)  GO  TO  Si 

A  VR 

3  1  bd 

Gl*UC**l5./b,J  *(3.ab9A('.8.15«9«3.e%*<U<f««a.-H,«0a3****3,)-l.:lb3»4F»»*/<» 

3  1  b  1 

*(5./b.) 

AV» 

3lba 

tFCul.Ol.o.J  , »wt , 

AvR 

31bo 

GO  1u  S3 

AvR 

3  1  b  4 

51 

Gl»«.d 

AvR 

3  1  b  5 

Go  II]  Si 

A  vR 

3  15b 

ba 

Gl*UI‘  •  •'S,/6.)*U<.259»M.e3S*«»0,Pd  4*0 **a, -0,  ^B«J*M»*3,J«*0**tS./b 

•  A  vR 

3  1  b  7 

•)*(l.»F|,FJJ3*U”O.<,3<4T»O«»£,*0tOO<1S»C»*3.) 

AVR 

31  bo 

IF"  (l»l  ,Lt.G,d)  Ut*d,3 

AVR 

3164 

53 

3X*BX*Gl*C4d*il"*'0(i.lJ 

AVR 

31  TO 

OS»OS*Om*i«SIL11 

AVR 

31  ’  1 

IF  (50,6.3,3)  GO  TO  S7 

AVR 

3i  Ta 

I  F  (*2  ,£u, 1 , ago . .o£ ,*« >  GO  TO  ST 

AVR 

3173 

V4»vH-u,5b*CN0*oT*«(5,/b,)*i)Fi««3(m) 

AVR 

•U7  4 

0o»OT/gk*5G£ 

AvR 

3  l  75 

«  *«F /oo 

AVR 

3  l  7  b 

IF (00,07,40)  uo  To  54 

AvR 

3177 

I  F  (k  ,  G£  ,  h  ,  5  ,  a  V) ,  J  ,  G£  ,  b  ,  5 )  Go  in  5* 

AvR 

3)7# 

IF (4. 1.7. 4, 5. *M). 0.1,7. 6. 5)  GO  TO  5b 

A  VR 

3179 

u8*u( •  •  (.5  .  /o  . )  *(d.aS9*n.B,,3*l»*0.Jld9*»»*a.-'3.rd4  3*rf**J.)*l,!'bS*0F. 

«  A  V  R 

3183 

• (b./b.) 

AvR 

3181 

IF (u0.lt,.).)  40*3. 

AvR 

3isa 

So  to  bo 

AVR 

3  10  3 

54 

G0*3.  <J 

A  vR 

3l8- 

GU  t!'  bo 

AVR 

«1  85 

55 

5i*o(''«A(b./o.)»tD.?S9*''.8k>5»9»;’,F!ob*R*Aa,*0.ndajAx*A3.)-«O«*l3,/b 

•  A  VR 

0  1  8b 

24 


A  V* 

3187 

I  F  (G2  •  It  •  o  •  •>  J  f-«f***,<* 

A  vR 

3188 

*56 

3  x  »  b  » -.Ga  *  t  o  2  •  u  "  »  «  a  ( L  1 ) 

A  V» 

0189 

0T*UT»i)h««iS  C L  1  ) 

A  VR 

3t9& 

57 

CUM  tNut 

A  V» 

8181 

m,*«f»*NGE«0i 

A  VS 

8192 

IF  (t/G.Lt. a.)  uS*r;S-u*«>«S (L 1 J 

A  V  R 

3l<»i 

Vi.«Vfl«vi,5s«CNa«(x*>*NGE-uS)*«(5a/0a)«(*“-Xi«2(H>»1200.),*3lLl5 

A  vfl 

8198 

3x»BX*Gl  »C92«  (*«i»iw2(Ll )  «  1U30,  )  «««S  (LI  ) 

A  VR 

8195 

os«os*(«w-itxaiLn  • 1000.0) 

AVR 

8198 

IF(K2,EU.O)  GO  TO  58 

a  ;» 

3197 

IF  (*2.£Q,  1  ,*nL>.wh0,l£.*><)  GU  TO  58 

AVR 

3198 

vxi»H»«,58»CM2*Ut*»(5,/8.)  (LI )  M000.)  »x3  (LI ) 

A  V  R 

3199 

bx*oX*G2*LN2*(HM»Xte2(Ll)*l880.)*"S(Ll] 

AVR 

8283 

0r»0T»(n»-i»2(Ll) »l J00.0) <>S(L1) 

A  vR 

0201 

56 

CUNTlNIJt 

AVR 

0282 

IF (*2.EU.0.*NU.L1 ,L£. J2)  GO  TU  81 

AVR 

3203 

IK Ll.LE.J»*IN,*N0,*a,t0.1)  GO  TO  81 

AVR 

0284 

60 

CUNHNUE 

AVR 

3235 

61 

Cunt iNut 

AVR 

0208 

IFUO.LT. 0.001)  GO  TO  <10 

AVR 

8207 

3x  »BX  «2 , l 1 7  «wx • «2, 

AVR 

320  3 

dI»tXP (4..HX1-1  , 

AVR 

8239 

Pv*«OvB»t»I«lACOS(7)«v«(l.-r««2.)««('.5)«r»0,01 

AVR 

82  U 

62 

C J*TINUt 

AVR 

8211 

?V0«PV0«lb,/PI 

AVR 

0212 

4* 

»!)•*#•»» 

AVR 

aaii 

FNl»l,-tXF(-i>,  •»«»•«,  5) 

AVR 

3214 

IF  (*O.Lf  .0.001)  GO  TO  <11 

AVR 

3215 

»V  IfcXP (4, aVK) -1 a 

AVR 

3218 

Gu»pvh/»v 

AVR 

8217 

v o  I  ■ v 0 I .GO ••0.5 

AVR 

3218 

41 

CUNTINUE 

AVR 

8219 

HtTuRN 

AVR 

3228 

2*0 

AVR 

3221 

I 
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